use "Data\MFPS Child Health Paper Replication Data - Child Level.dta", clear

* Set Panel Variables
xtset caseid_n year

*Set Scheme
set scheme modern


* Create global of covariates 
global covariates "i.cluster_num work_bl ever_use_bl curr_use_bl sex_age_bl i.age_group_bl edu_primary_bl i.birth_order religion_r_bl ethnicity_r_bl i.total_birth_bl child_sex"

*Create global of age covariates
mkspline age1 6 age2 12 age3 18 age4 24 age5 30 age6 = child_age_months
global agevars "i.chi1m age1-age6"

forvalues i = 1/6{
	gen chi_age_group_`i' = age`i'>0
}

*Covariates that can only be used with the first child
*birthweight timetobreast wantedness 
* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* %%%%% REGRESSIONS AND TABLES %%%%%
* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

/////////////////////////////////////
** %%% CONTROL MEANS %%% **
/////////////////////////////////////

tabstat haz06 stunted stunted_ex whz06 child_hem if year==2017 & treatment==0 & index_dummy==1 & child_age_months>6, statistic(mean) columns(statistics)

tabstat credi_score credi_z_score if year==2018 & treatment==0 & index_dummy==1 & child_age_months>6, statistic(mean) columns(statistics)

tabstat haz06 stunted stunted_ex whz06 waz06 child_hem if year==2017 & index_dummy == 1 & born_during_study==0, statistic(mean) columns(statistics)

tabstat credi_score credi_z_score if year==2018 & index_dummy == 1 & born_during_study==0, statistic(mean) columns(statistics)

/////////////////////////////////////
** %%% OLS ESTIMATION %%% **
/////////////////////////////////////
*Store P values for MHT
putexcel set "Results\ITT\\MHT p-values.xlsx", replace
	putexcel A1= ("Variable")
	putexcel B1= ("All Index unadjusted")
	putexcel C1 = ("All Index age adjusted")
	putexcel D1 = ("All Index full adjusted")
	putexcel E1 = ("Pre-index unadjusted")
	putexcel F1 = ("Pre-index age adjusted")
	putexcel G1 = ("Pre-index full adjusted")

*Main results local
local outcomes haz06 stunted stunted_ex whz06 waz06 child_hem credi_score credi_z_score
local i = 1
local cell = 2
foreach y of local outcomes{
	if inrange(`i',1,6) {
		local year = 2017
	}
	if inrange(`i',7,8){
		local year = 2018
	}
	local name: variable label `y'
	putexcel A`cell' = ("`name'")
	estimates clear

	*Without children in the exclusive breastfeeding period
	eststo a1: qui reg `y' treatment if year==`year' & index_dummy==1 & ///
		child_age_months>6, cluster(woman_id)
			local p = r(table)[4,1]
			putexcel B`cell' = (`p')
			
	eststo a2: qui reg `y' treatment $agevars if year==`year' & index_dummy==1 & ///
		child_age_months>6, cluster(woman_id)
			local p = r(table)[4,1]
			putexcel C`cell' = (`p')
			
	eststo a3: qui reg `y' treatment $agevars $covariates if year==`year' & ///
		index_dummy==1 & child_age_months>6, cluster(woman_id)
			local p = r(table)[4,1]
			putexcel D`cell' = (`p')
			
	eststo a4: qui reg `y' treatment if year==`year' & index_dummy==1 ///
		& born_during_study==0, cluster(woman_id)
			local p = r(table)[4,1]
			putexcel E`cell' = (`p')
			
	eststo a5: qui reg `y' treatment $agevars if year==`year' & index_dummy==1 ///
		& born_during_study==0, cluster(woman_id)
			local p = r(table)[4,1]
			putexcel F`cell' = (`p')
			
	eststo a6: quietly reg `y' treatment $agevars $covariates if year==`year' & ///
		index_dummy==1 & born_during_study==0, cluster(woman_id)
			local p = r(table)[4,1]
			putexcel G`cell' = (`p')

	*** Tabulate Results
	esttab a1 a2 a3 a4 a5 a6 ///
	using "Results\ITT\\`name' (ITT).rtf", ///
	replace se ///
	scalar(N r2 F)  mlabels("All Index" "All Index" ///
	"All Index" "Index & Born Prior" "Index & Born Prior" "Index & Born Prior") collabels(none) ///
	starlevels(* .1 ** .05 *** .01)     ///
	compress keep(treatment _cons) ///
	order(treatment)	

	local i = `++i'
	local cell = `++cell'
}


/////////////////////////////////////////
** %%% Inverse Probability Weight %%% **
/////////////////////////////////////////

local outcomes haz06 stunted stunted_ex whz06 waz06 child_hem credi_score credi_z_score
local i = 1
foreach y of local outcomes{
	if inrange(`i',1,6) {
		local year = 2017
		local lag = "" 
	}
	if inrange(`i',7,8){
		local year = 2018
		local lag = "2"
	}
	local name: variable label `y'
	
	*Generare measurement indicators
	gen measured_`year' = (`y'!=.) if year == `year'
	replace measured_`year' = F`lag'.measured_`year' if year == 2016
	replace measured_`year' = 0 if year == 2016 & measured_`year'==.
	gen measured_2016 = (`y'!=.) if year == 2016
	replace measured_2016 = L`lag'.measured_2016 if year== `year'
	
	
	*Estimate propensity scores 
	eststo pw_`y': logit measured_`year' treatment measured_2016 $agevars $covariates if year == 2016 & born_during_study == 0 & index_dummy==1, robust or
	predict measured_reg if year ==2016
	replace measured_reg = L`lag'.measured_reg if year == `year'
	
	*Plot Propensity Scores
	twoway (kdensity measured_reg if year == 2016 & index_dummy==1 & born_during_study==0 & measured_`year' == 0, lpattern(dash) lcolor(gray)) ///
	(kdensity measured_reg if year == 2016 & index_dummy==1 & born_during_study==0 & measured_`year' == 1, lcolor(maroon)), ///
	title(`name') xtitle("Pr(Measurement at First-Year Follow-Up (FYP))") ytitle("Density") legend(order(1 "Not Measured at FYP" 2 "Measured at FYP") pos(1) ring(0))
	graph export "Results\Figures\Propensity Scores\\`name'.png", replace
	
	*Estimate weighted least squares
	eststo wls_`y': qui reg `y' treatment $agevars $covariates if year==`year' & index_dummy==1 ///
		& born_during_study==0 [aw = 1/measured_reg], cluster(woman_id)
		
local i = `++i'
cap drop measured_`year' measured_2016 measured_reg
		
}

esttab pw_* using "Results\Propensity Scores\First Stage.csv", ///
replace se ///
scalar(N r2 F) collabels(none) ///
starlevels(* .1 ** .05 *** .01)     ///
compress eform
	
esttab wls_* using "Results\Propensity Scores\Second Stage.csv", ///
replace se ///
scalar(N r2 F) collabels(none) ///
starlevels(* .1 ** .05 *** .01)     ///
compress 

//////////////////////////////////////////
** %%% Heckman Selection %%% ** 
//////////////////////////////////////////
local outcomes haz06 stunted stunted_ex whz06 waz06 child_hem credi_score credi_z_score
local i = 1
foreach y of local outcomes{
	if inrange(`i',1,6) {
		local year = 2017
	}
	if inrange(`i',7,8){
		local year = 2018
	}
	local name: variable label `y'
	
	eststo h_`y': heckman `y' treatment $agevars $covariates if year==`year' & born_during_study==0 & index_dummy==1, select(i.w1_case_enumid $agevars $covariates)
	
	*joint test of selection FE 
	gen select = 0 
	replace select = 1 if `y' != . & year==`year' & born_during_study==0 & index_dummy==1
	
	qui reg select i.w1_case_enumid $agevars $covariates if year==`year' & born_during_study==0 & index_dummy==1
	testparm i.w1_case_enumid
	drop select
	
	local i = `++i'
	
}

esttab h_* using "Results\Heckman Selection\heckman selection results.csv", ///
replace se ///
scalar(N r2 F) collabels(none) ///
starlevels(* .1 ** .05 *** .01)     ///
compress 

		
/////////////////////////////////////////
** %%% Kling-Leibman Bounds %%% **
/////////////////////////////////////////
preserve
keep if born_during_study==0 & index_dummy==1
tsset caseid_n year
tsfill, full

global fill "treatment chi1m cluster_num work_bl ever_use_bl curr_use_bl sex_age_bl age_group_bl edu_primary_bl birth_order religion_r_bl ethnicity_r_bl total_birth_bl child_sex woman_id index_dummy born_during_study"
foreach i of global fill{ 
    replace `i' = L.`i' if year ==2017 & `i' == . 
	replace `i' = L2.`i' if year ==2018 & `i' == . 
}

replace child_age_months = L.child_age_months+12 if year ==2017 & child_age_months==.
replace child_age_months = L2.child_age_months+24 if year ==2018 & child_age_months==.
drop age1-age6
mkspline age1 6 age2 12 age3 18 age4 24 age5 30 age6 = child_age_months


local outcomes haz06 stunted stunted_ex whz06 waz06 child_hem credi_score credi_z_score
local i = 1
foreach y of local outcomes{
	if inrange(`i',1,6) {
		local year = 2017
	}
	if inrange(`i',7,8){
		local year = 2018
	}
	local name: variable label `y'
	
	*Generate means and standard deviations
	summ `y' if year == `year' & born_during_study==0 & index_dummy==1 & treatment==0
	local mu_c = r(mean)
	local stand_c = r(sd)
	
	summ `y' if year == `year' & born_during_study==0 & index_dummy==1 & treatment==1
	local mu_t = r(mean)
	local stand_t = r(sd)
	
	**Gen bounding outcomes
	gen low20 = `y'
	replace low20 = `mu_c'+0.2*`stand_c' if low20 == . & treatment == 0 
	replace low20 = `mu_t'-0.2*`stand_t' if low20 == . & treatment == 1 
	
	gen low10 = `y'
	replace low10 = `mu_c'+0.1*`stand_c' if low10 == . & treatment == 0 
	replace low10 = `mu_t'-0.1*`stand_t' if low10 == . & treatment == 1 
	
	gen high20 = `y'
	replace high20 = `mu_c'-0.2*`stand_c' if high20 == . & treatment == 0 
	replace high20 = `mu_t'+0.2*`stand_t' if high20 == . & treatment == 1 
	
	gen high10 = `y'
	replace high10 = `mu_c'-0.1*`stand_c' if high10 == . & treatment == 0 
	replace high10 = `mu_t'+0.1*`stand_t' if high10 == . & treatment == 1 
	
	*Conduct Regressions
	eststo low20: quietly reg low20 treatment $agevars $covariates if year==`year' & ///
		index_dummy==1 & born_during_study==0, cluster(woman_id)
		
	eststo low10: quietly reg low10 treatment $agevars $covariates if year==`year' & ///
		index_dummy==1 & born_during_study==0, cluster(woman_id)
		
	eststo high10: quietly reg high10 treatment $agevars $covariates if year==`year' & ///
		index_dummy==1 & born_during_study==0, cluster(woman_id)
		
	eststo high20: quietly reg high20 treatment $agevars $covariates if year==`year' & ///
		index_dummy==1 & born_during_study==0, cluster(woman_id)
		
	*** Tabulate Results
	esttab low20 low10 high10 high20 ///
	using "Results\ITT\Bounds\\`name' (bounds).csv", ///
	replace se ///
	scalar(N r2 F)  mlabels("Low 20" "Low 10" "High 10" "High 20") collabels(none) ///
	starlevels(* .1 ** .05 *** .01)     ///
	compress keep(treatment _cons) ///
	order(treatment)
	
	local i = `++i'
	cap drop low20 low10 high10 high20
	estimates clear

}
restore

//////////////////////////////////////////////////////
* CREATE FIGURES 
/////////////////////////////////////////////////////
twoway (kdensity haz06 if year ==2016 & treatment==0 & index_dummy==1 & born_during_study==0, lpattern(dash) lcolor(gray)) ///
(kdensity haz06 if year ==2016 & treatment==1 & index_dummy==1 & born_during_study==0, lcolor(maroon)), ///
title("2016") ytitle("Density") xtitle("Height-for-Age Z-Scores") legend(order(1 "Control" 2 "Treatment")) 
graph export "Results\Figures\HAZ Density 2016.png", replace

twoway (kdensity haz06 if year ==2017 & treatment==0 & index_dummy==1 & born_during_study==0, lpattern(dash) lcolor(gray)) ///
(kdensity haz06 if year ==2017 & treatment==1 & index_dummy==1 & born_during_study==0, lcolor(maroon)), ///
title("2017") ytitle("Density") xtitle("Height-for-Age Z-Scores") legend(order(1 "Control" 2 "Treatment")) 
graph export "Results\Figures\HAZ Density 2017.png", replace

twoway (kdensity credi_z_score if year ==2018 & treatment==0 & index_dummy==1 & born_during_study==0, lpattern(dash) lcolor(gray)) ///
(kdensity credi_z_score if year ==2018 & treatment==1 & index_dummy==1 & born_during_study==0, lcolor(maroon)), ///
ytitle("Density") xtitle("CREDI Z-Scores") legend(order(1 "Control" 2 "Treatment")) 
graph export "Results\Figures\CREDI Density.png", replace

//////////////////////////////////////////////////////
* MHT Corrections
/////////////////////////////////////////////////////

** Multiple Hypothesis Testing 
import excel "Results\ITT\\MHT p-values.xlsx", sheet("Sheet1") firstrow clear

*Create Sharped Q-Values 
qqvalue AllIndexunadjusted, method(simes) qvalue(Q_AllIndexunadjusted)
qqvalue AllIndexageadjusted, method(simes) qvalue(Q_AllIndexageadjusted)
qqvalue AllIndexfulladjusted, method(simes) qvalue(Q_AllIndexfulladjusted)
qqvalue Preindexunadjusted, method(simes) qvalue(Q_Preindexunadjusted)
qqvalue Preindexageadjusted, method(simes) qvalue(Q_Preindexageadjusted)
qqvalue Preindexfulladjusted, method(simes) qvalue(Q_Preindexfulladjusted)
	
export excel "Results\ITT\\MHT p-values.xlsx", sheet("Sheet1") firstrow(variables) replace 
